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While restricted single-reference coupled cluster theory truncated to singles and doubles (CCSD) 
provides very accurate results for weakly correlated systems, it usually fails in the presence of 
static or strong correlation. This failure is generally attributed to the qualitative breakdown of 
the reference, and can accordingly be corrected by using a multi-determinant reference, including 
higher-body cluster operators in the ansatz, or allowing symmetry breaking in the reference. None 
of these solutions are ideal; multi-reference coupled cluster is not black box, including higher-body 
cluster operators is computationally demanding, and allowing symmetry breaking leads to the loss of 
good quantum numbers. It has long been recognized that quasidegeneracies can instead be treated 
by modifying the coupled cluster ansatz. The recently introduced pair coupled cluster doubles 
(pCCD) approach is one such example which avoids catastrophic failures and accurately models 
strong correlations in a symmetry-adapted framework. Here we generalize pCCD to a singlet- 
paired coupled cluster model (CCDO) intermediate between coupled cluster doubles and pCCD, 
yielding a method that possesses the invariances of the former and much of the stability of the 
latter. Moreover, CCDO retains the full structure of coupled cluster theory, including a fermionic 
wave function, antisymmetric cluster amplitudes, and well-defined response equations and density 
matrices. 


I. INTRODUCTION 

Over the past few decades, single-reference coupled 
cluster (CC) theorjff® has reigned unchallenged as 
the dominant paradigm for the accurate description of 
weakly correlated systems in quantum chemistry. This 
dominance is justly earned. Coupled cluster delivers very 
accurate results both for energies and for other proper¬ 
ties, with reasonable polynomial scaling. With an ade¬ 
quate reference determinant, it is size consistent (correct 
separability) and size extensive (correct scaling with sys¬ 
tem size). The method is capable of treating closed shells 
and open shells, ground states and excited states. 

Where standard forms of single-reference coupled clus¬ 
ter fail is in the description of problems where strong or 
static correlation effects are important. In such cases, 
the failure of single-reference coupled cluster theory is 
almost always attributed to the reference determinant: 
the presence of static correlation, we are told, means that 
the restricted Hartree-Fock (RHF) reference is not even 
qualitatively reasonable. This is, of course, true. But 
why should the exponential coupled cluster wave func¬ 
tion fail catastrophically, as we see, for example, in the 
dissociation of N 2 ,® or in the large U limit of the Hubbard 
Hamiltonian,® or in the large G limit of the attractive 
pairing Hamiltonian ® 7 To be sure, a correlated wave 
function built atop a poorly chosen reference need not 
give accurate predictions, but we would like the predicted 


energy to at least be real and finite! In other words, one 
must tame the failures of coupled cluster theory built 
atop a symmetry adapted reference. 

One way to overcome these failures is simply to work 
harder. That is, one can add higher and higher lev¬ 
els of excitation to the cluster operator or at least se¬ 
lected higher excitations; eventually, the cluster operator 
will contain excitations of sufficiently high rank to over¬ 
come the deficiencies of the reference determinant. After 
all, coupled cluster theory is formally exact when one 
includes all possible excitations in the cluster operator. 
Multireference methods, in other words, are not needed 
for strongly correlated systems, if one can use a suffi¬ 
ciently complete single-reference method. Put this way, 
the failures of single reference coupled cluster theory for 
strongly correlated systems are failures of computational 
ability, not deficiencies of the theory. 

The purpose of this manuscript is to emphasize that 
coupled cluster theory can be protected from these catas¬ 
trophic failures by working smarter instead. We seek im¬ 
provement, that is, by removal rather than by addition. 
This is not a new idea. Numerous authors have proposed 
methods that yield reasonable predictions in the presence 
of near degeneracies by modifying coupled cluster theory 
or variants thereof.®®! We would like a simple concep¬ 
tual framework which allows us to identify the offending 
terms in traditional coupled cluster theory so that by ex¬ 
cluding them we can avoid the downfall of the method. 
Of course, this is easier said than done, but we here show 
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one particular case which substantially removes the typ¬ 
ical failures in repulsive Hamiltonians. 

What causes these failures? Closure of the occupied- 
virtual gap is commonly cited as the culprit, but this 
cannot be the whole story. It is, of course, true that 
closure of the gap causes perturbation theory to fail, but 
the infinite order summations in coupled cluster theory 
can overcome a vanishing gap. For example, the coupled 
cluster doubles (CCD) equations can schematically be 
written as 


Ae-t = v + f(t) (1) 

where v is a two-electron integral; closure of the gap 
(Ac —► 0) just means that we must have t such that 
v + f(t) = 0. Thus, for example, CCD is well-behaved 
for the dissociation of H 2 despite the gap closing; simi¬ 
larly, we can obtain perfectly reasonable results for the 
homogeneous electron gas even though it is metallic, 21 
or for the repulsive pairing Hamiltonian even when all 
occupied levels have higher energy than all virtual levels. 

If the problem is not the vanishing gap, what else 
could it be? The key insight can be obtained by con¬ 
sidering the underlying mean-field theory. As correla¬ 
tions become stronger, RHF becomes unstable toward a 
symmetry broken unrestricted Hartree-Fock (UHF). The 
symmetry spontaneously (but artificially) breaks at the 
mean-field level, and one can safely build coupled cluster 
wave functions starting from the UHF reference at the 
cost of good quantum numbers (symmetry labels) which 
are difficult to recover once lost. 

It seems clear that the failures of single-reference cou¬ 
pled cluster theory are in some way tied to these in¬ 
stabilities. After all, the particle-hole random phase 
approximation (RPA) lurks within the coupled cluster 
equationsp2H23 an d the RPA delivers unphysically large 
correlation energies as RHF approaches the instability 
point and complex correlation energies past it. The prob¬ 
lems caused by these instabilities are at least partially 
screened even at the simple CCD level, as is clear by 
the fact that CCD can dissociate H 2 correctly where 
RPA cannot. However, fully eliminating the problems 
caused by these instabilities would appear to require go¬ 
ing beyond CCD in general. This we are loathe to do, 
simply because it is computationally demanding. As 
we have said, we would prefer an alternative frame¬ 
work which does not greatly increase computational com¬ 
plexity. Much of the inspiration in the present paper 
comes from the recently introduced pair CCD (pCCD) 
metliod j 20 l 25 lESI] which seems immune to problems with 
instabilities and provides a reasonable description of 
strong correlation effects. However, pCCD has its own 
limitations and requires a full orbital optimization, which 
is far from being black-box. Our purpose here is to elim¬ 
inate as many of the vices of pCCD as possible while 
retaining most of its virtues. 

In the remainder of this paper, we will first illustrate 
how the failures of CCD can be avoided by removing the 
terms sensitive to instabilities at the RPA level. This 


approach, however, is not ideal, because leaving out en¬ 
tire classes of terms removes the association between the 
cluster operator on the one hand and an approximate 
solution to the Schrodinger equation on the other. We 
remedy this deficiency in Sec. Ill by removing parts of 
the cluster operator itself, thereby retaining the eigenvec¬ 
tor property, and show several illustrative results before 
offering concluding remarks. 


II. FAILURES IN COUPLED CLUSTER 
DOUBLES 


In this section we wish to provide a few pedagogical ex¬ 
amples which illustrate clearly how failures of restricted 
CCD can be related to symmetry instabilities in the un¬ 
derlying reference determinant and can be ameliorated 
by simply deleting terms from the amplitude equations. 
To do so, we should say a few words about the Hartree- 
Fock stability conditions. 

Solving the Hartree-Fock equations guarantees only 
that the energy is stationary with respect to orbital rota¬ 
tions. To obtain a local minimum, we also need that the 
Hartree-Fock orbital Hessian has only positive eigenval¬ 
ues. If the Hessian has a negative eigenvalue, the Hartree- 
Fock solution is unstable and a lower energy solution can 
be found. Particle-hole RPA (ph-RPA) is closely related 
to the particle-hole Hartree-Fock Hessian, and when this 
Hessian has a negative eigenvalue, ph-RPA predicts a 
complex correlation energy. Similarly, particle-particle 
(pp-RPA) is related to a Hartree-Fock Hessian which is 
sensitive toward a number-broken mean-field reference; 
when the Hartree-Fock is unstable to number symme¬ 
try breaking, the Hessian develops a negative eigenvalue 
and the pp-RPA correlation energy becomes complex. 
One can of course use an eigenvector corresponding to 
a negative eigenvalue to find a mean-field solution with 
lower energy and typically lower symmetry. If the mean- 
field breaks symmetry, the negative eigenvalue disap¬ 
pears in favor of a Goldstone mode with zero eigenvalue 
in the broken-symmetry Hessian; thus, the RPA built 
atop a symmetry-broken reference yields a real correla¬ 
tion energyPH We should emphasize that the RPA we 
discuss here includes the full exchange interaction. By 
neglecting exchange to form direct RPA, one severs the 
link between RPA and the Hartree-fock stabil ity pr oblem 
and thereby obtains real correlation energies!— I— I 

While CCD seems a world removed from stability 
of Hartree-Fock solutions, contact is made through the 
RPA. Both ph-RPA and pp-RPA can be cast as chan¬ 
nels of CCD, 22 ^ so it seems at least plausible that the 
Hartree-Fock instabilities which cause the failures of RPA 
may be ultimately responsible for the failures of CCD 
as well. This in turn suggests that these failures might 
be remedied by deleting the appropriate RPA-like terms 
from the CCD amplitude equations. 

Let us take a moment to be precise. In CCD, the wave 
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function is given as 

l*>=e Ta |0> (2) 

where |0) is a reference determinant (taken here to always 
be the RHF determinant) and the cluster operator T 2 
creates double excitations: 

T 2 = \j2 c “ c b c o c i ( 3 ) 

ijab 

where spinorbitals indexed by i, j. k , l (a, b , c, d) are 
occupied (unoccupied) in |0). The correlation energy is 

e = < 4 > 

ijab 

where vj b = (ij\\ab) is an antisymmetrized two-electron 
integral in Dirac notation, while the wave function am¬ 
plitudes are obtained by solving 

n _ fCL ±cb . fb +ac rk j.ab rk ±ab , ab /r\ 

^ Jc ''ij ' Jc ^ij Ji ^kj Jj ''ik ' ^ij v'-v 

_L_ 1 kl , 1 j cd ab , \ f cd + ab kl 

' 2 “ij ' 2 ^ij ^cd ' ^ ''ij ^kl “cd 

i >ac kb , >cb n .ka , ±ac j.bd kl 

' ^ik V c j ' ^kj vci ' T'ik ^jl ^cd 

j.bc n .ka j.ca kb j.bc +ad kl 
^ik u cj ^kj Vci ^ik c jl ^cd 

^ kl (a. cd j.ab | ±cd j.ab , ±ad j.cb , ±db j.ac\ 

2 ^cd \ l il ^kj ' ''Ij ''ik ' ^kl ^ij ' ^kl ^ij ) 

in terms of the Fock operator = ( q\f\p ); in this equa¬ 
tion, summation over repeated indices is implied. We 
refer to the terms on the first line of the CCD ampli¬ 
tude equation as the driver; the terms on the next three 
lines we call the ladder, ring, and crossed-ring terms, and 
we call the terms on the last line mosaic termsP^ One 
can obtain ph-RPA from CCD simply by retaining only 
the driver and ring terms (modulo an overall factor of 
two on the correlation energy which in this work we will 
neglect), and pp-RPA by retaining only the driver and 
ladder terms. 

All of this means that if the RPA instabilities drive 
failures of CCD, one should obtain well-behaved results 
by eliminating the ladder terms in systems for which pp- 
RPA exhibits instabilities and by eliminating the ring 
terms in systems for which ph-RPA becomes unstable; in 
the latter case, one should also remove the crossed-ring 
terms so as to retain proper fermionic antisymmetry of 
the amplitudes f“ b . To see that eliminating the offend¬ 
ing RPA terms does indeed cure the failures of CCD, we 
provide several illustrative examples. 

Consider, then, the dissociation of the N 2 molecule. 
We will work in the STO-3G basisp^ which is useless for 
coupled cluster calculations in general but which here 
serves to isolate the effects of static correlation fairly 
well as we have essentially a valence-only calculation. 
As Fig. [T] shows, CCD breaks down very badly in this 
case, predicting a dissociation limit very close to the equi¬ 
librium energy. The ph-RPA energy is disastrous ev¬ 
erywhere, and becomes complex near equilibrium. On 



FIG. 1. Total energies of the N 2 molecule in the STO-3G basis 
set as a function of bond length. Past the point shown, ph- 
RPA yields a complex correlation energy. “lm-CCD” denotes 
CCD where we have eliminated ring and crossed-ring terms. 



FIG. 2. Total energies in the 6-site Hubbard Hamiltonian 
at half filling as a function of interaction strength U /t. The 
ph-RPA yields a complex correlation energy past the point 
shown. “lm-CCD” denotes CCD where ring and crossed-ring 
terms have been excluded. 

the other hand, removing ring and crossed-ring contrac¬ 
tions to obtain what we have called lmCCE^ produces a 
well-behaved dissociation curve (in that the energy rises 
monotonically from equilibrium to dissociation), albeit 
one which is by no means energetically satisfactory. 

The same happens in the repulsive Hubbard Hamilto¬ 
nian, which is given by 

H = -tJ2 c L c ^+ u n Pt n Pi ( 6 ) 

(pq),a P 

where p and q index lattice sites and a indexes spins, 
n p<T = Cp^ c Pa is a spinorbital number operator, and the 
notation ( pq) indicates that sites p and q are nearest 
neighbors; we take the lattice to be periodic so that sites 
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FIG. 3. Total energies in the attractive pairing Hamiltonian 
as a function of interaction strength G. Past the points shown, 
pp-RPA and CCD yield a complex correlation energy, “rxm- 
CCD” denotes CCD where laddter terms have been omitted. 


1 and N in an iV-site Hubbard Hamiltonian are adjacent. 
In Fig. [2j we see that CCD fails again for the 6-site Hub¬ 
bard model at half-filling. While the exact energy rises 
monotonically to zero as U approaches infinity, the CCD 
energy quickly turns over and gets increasingly negative 
as U gets larger, and ph-RPA yields a complex correlation 
energy even for relatively small U. Again, removing the 
ring and crossed-ring terms cures this qualitative failure 
without, however, providing reasonable energies. 

Lastly, consider the attractive pairing Hamiltonian, 
given by 

H = € p n v„ - G 4t c k c n c ?t > ( 7 ) 

p,cr pq 

where G > 0 and the levels are equally spaced (e p = p). 
For this problem, there is an instability toward number 
symmetry breaking. Both pp-RPA and in fact CCD even¬ 
tually predict complex correlation energies, but as can be 
seen from Fig. [3j removing the ladder terms in what we 
have denoted rxmCCD provides a stable but inaccurate 
curve. 

It is clear that we can indeed cure the failures of CCD 
by simply deleting entire classes of terms from the ampli¬ 
tude equations. The cure, however, is often worse than 
the disease. We shall do better shortly, but it is instruc¬ 
tive to see why eliminating terms from the CCD equa¬ 
tions may yield poor results. 

The CCD energy and amplitude equations shown ear¬ 
lier are given by 

-Eccd = (0|e -T2 H e T2 |0), (8a) 

0= ($“ h | e - T2 He T2 |0), (8b) 

where the state <!>“*) is a doubly-excited determinant. It 
is clear that CCD can be viewed as writing a similarity- 
transformed Hamiltonian H = exp(—T 2 ) H exp(T 2 ) such 


that the reference determinant |0) is a right-hand eigen¬ 
vector within the space of the reference determinant and 
all double excitations. In other words, the CCD ampli¬ 
tude equations are obtained by solving the Schrodinger 
equation for a transformed Hamiltonian within a sub¬ 
space. This is no longer the case when one removes terms 
whole cloth. In other words, when one eliminates entire 
classes of terms from the CCD amplitude equations, one 
is no longer solving the Schrodinger equation; accord¬ 
ingly, it should not be too surprising that the results 
thereby obtained are unsatisfactory. 


III. SINGLET-PAIRED CCD 

Ideally, we would like to formulate CCD in terms of 
a cluster operator which is insensitive to instabilities in 
the reference without having to remove entire classes of 
terms, because by doing so, we guarantee that we have 
a solution to the Schrodinger equation at least within a 
subspace of the full Hilbert space of the problem. While 
it is not entirely clear how one might accomplish this 
in a general way, we do have a specific example of a 
CCD model which avoids these problems. That model 
is pCCD, in which the cluster operator is rewritten as 
simply 

Tf GB = J2 T i c U c l c h c n ( 9 ) 

ia 

where here and for the remainder of this manuscript or¬ 
bital labels i, j, a, b, and so forth refer to spatial or¬ 
bitals. Thus, for example, c' a cj Ci x c,; t removes two elec¬ 
trons from spatial orbital i and places them in spatial 
orbital a. Remarkably, this simplification leads to a well- 
behaved model which neither overcorrelates wildly nor 
predicts a complex correlation energy in systems which 
are unstable toward a UHF reference. To some extent, 
we can attribute the success of pCCD in this regard to 
the fact that it can include the perfect pairing valence 
bond and antisymmetrized product of strongly orthogo¬ 
nal geminals wave functions, albeit with coefficients de¬ 
termined by similarity transformation rather than via the 
varitional principle. 20 

However, while restricted pCCD provides a reasonable 
accounting for the paired correlations, it does not include 
any correlations between electrons in different spatial or¬ 
bitals. Moreover, the result of a pCCD calculation is not 
invariant to rotations amongst the occupied orbitals or 
amongst the virtual orbitals, which necessitates an or¬ 
bital optimization which can be somewhat tedious. The 
question then becomes whether we can eliminate these 
flaws without reintroducing the failures for strongly cor¬ 
related systems. 

In order to do so, however, it will prove useful to first 
revisit standard CCD in the RHF framework, for which 
the cluster operator is singlet spin adapted and can be 
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FIG. 4. Left panel: Total energies of the N 2 molecule in the STO-3G basis set as a function of bond length. Right-panel: Total 
energies in the 6-site Hubbard Hamiltonian at half filling as a function of interaction strength UIt. In both cases, “lm-CCD” 
denotes CCD where particle-hole contractions have been excluded. 


conveniently written a: 


E6i 


1 


T * = *'Z I u'E 


f f 

Cl cl Cj c Ci 

a a Of J£, L o 


( 10 ) 


ijab 




It is important to note that as a two-body operator, T 2 
actually contains two independent singlet components 
which we have parameterized by the single set of am¬ 
plitudes T“ b = T ba = f°t 6 ij ii | 36 | 37 | That i s we have 

A L J J b 


t 2 = tJ 0] + t] 1] , 


Tr 


[ 0 ] 


= - V CT ob P f 
9 l 3 < 


ab J 


ijab 


T 


[ 1 ] 


= 9 


(lla) 

(lib) 

(11c) 


ijab 


where the pair operators P a b and Q a b are given by 


Pij 

— ^2 ( c Jt Ci i + Ci t c h) > 

(12a) 

Qij 

= Cj f Ci T , 

(12b) 

Qij 

= c ii c h 1 

(12c) 

Qij 

= \/i ^ C,7t Cii Cif 

(12d) 


and where the notation Q' ab ■ Qij is shorthand for 

Ql-Qij = (Qlb) f Qtj + iQab) 1 QTj + {Qab) 1 Q%- (is) 

Notice that Pij is symmetric on interchange of i with j 

t.Vnis r, . . ....... .. . . 

f-3 13 


whereas Qij is antisymmetric; thus, cr“ b and 7r“ b obey 


respectively 


tpab _i_ rp, 

ab _ ba _ ab _ ba _ ij ' ij 

a ij ~ a ij ~ a ji ~ a ji ~ 2 


ba 


rpab _ rp, 

_ab ba ab ba ij ij 

’! I • • — 1 1 • ■ — j i • • — y 1 • • — 

11 IJ IJ Jl Jl 2 


ba 


(14a) 

(14b) 


so that aij (7T? b ) is the part of T“ b symmetric (antisym¬ 
metric) on interchange of i with j. The operators P % j 
are singlet-paired, while Q ^ are triplet paired; in other 
words, the wave function created by pl|—) is a singlet, 
while those created by the three components of Qf|—) 
are triplets, where |—) is the physical vacuum. Note that 
t] 0 ^ correlates only electrons of opposite spin, while T ^ 
includes both opposite-spin and same-spin correlation. 
However, both T. ^ and yield contributions from all 
channels of CCD, i.e. they retain ladder, ring, crossed- 
ring, and mosaic terms. 

Clearly, pCCD eliminates the triplet-paired operators 
entirely and retains only the diagonal i = j and a = b 
portions of the singlet-paired operators. In other words, 
in pCCD we eliminate entirely and retain only the 

diagonal portion of T. ^. We could add correlations be¬ 
tween electrons in different spatial orbitals and restore 
invariance to occupied-occupied and virtual-virtual mix¬ 
ing by simply relaxing the restriction that we include only 
the diagonal portions of the singlet-paired operators. Put 
differently, we can generalize pCCD to what we will call 
CCDO by simply writing 

T 2 CCD0 = if 1 . (15) 

Operationally, this just means that in a standard RHF- 
based CCD code which uses t/T 1 as the fundamen- 
tal ingredient, as many codes do,® one needs only to 
symmetrize the amplitudes T“ b iteration by iteration; 
in other words, one makes the replacement T“ b —> 
1/2 (T“ b + 7^“) every iterative cycle. 

Figure [4] shows the results of CCDO for the dissociation 
of N 2 and for the Hubbard Hamiltonian. It is clear that 
CCDO remedies the failures of CCD while still including 
much of the correlations that are lost by eliminating the 
ring terms in CCD. These examples, however, are fairly 
trivial, so we will present several more results to establish 
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FIG. 5. Geometry of H 4 on a circle. 



FIG. 6 . Total energies of H 4 as a function of angle 0. 


the point and to show how one can build atop CCDO in 
a natural way. 

A strict test for coupled cluster is the case of four hy¬ 
drogen atoms symmetrically distributed on a circle of 
radius R = 1.738A, as depicted in Fig. [5p! For small or 
large angles 0 , the system resembles two H 2 units that 
are reasonably well separated, but as 0 passes through 
90°, the four atoms form a square and there is a degener¬ 
acy. The exact energy is smooth as a function of 0, but at 
the RHF level there is a cusp at 90°, just as with the ro¬ 
tation of the C-C bond in ethylene. Figure [ 6 ] sho ws total 
energies for H 4 in Dunning’s DZP basis set ® 41 all built 
atop the RHF reference. Both CCD and coupled clus¬ 
ter with singles and doubles (CCSD) have a cusp and 
incorrectly predict a local minimum at 90°. While the 
cusp remains in CCDO and in CCSDO ( i.e. CCDO with 
the full inclusion of single excitations), it is smoothed 
out somewhat and both CCDO and CCSDO correctly pre¬ 
dict a maximum in the energy at 90° instead. Qualita¬ 
tively, then, CCDO and CCSDO are much closer to correct 
than are CCD and CCSD; CCSDO is fortuitously semi- 
quantitative for this problem but not, as we shall see, in 



r n-n (A) 

FIG. 7. Total energies of the N 2 molecule in the cc-pVDZ 
basis set as a function of bond length. Where CCSD and 
CCSDT break down for large R, CCSDO and CCSDTO are 
properly stable. 


general. 

Next, we consider N 2 with the slightly more complete 
cc-pVDZ basis setas seen in Fig. [?] The failure of 
CCD (and therefore of CCSD) is lessened in the larger 
basis, which we can qualitatively understand as arising 
from the much larger number of virtual orbitals allowing 
for significantly more dynamic correlation. Nonetheless, 
CCSD breaks down in the larger basis where CCSDO does 
not. One could also employ Brueckner orbitals in a CCDO 
analogue of Brueckner doubles, which we would denote 
as BDO, but results do not differ much from those of 
CCSDO, so we have not shown them here. 

These calculations permit us also to make a second 
point. In traditional coupled cluster theory, only the Tf 
and X 2 equations are nonlinear in their respective opera¬ 
tors, i.e. the T 3 equations are linear in T 3 (though they 
contain terms like H T 2 J 3 ), the X 4 equations are linear in 
T 4 , and so on. If the failures in CCD can be attributed to 
the terms quadratic in T 2 , therefore, it seems reasonable 
to suppose that including higher-order cluster operators 
atop CCDO would immediately lead to a well-behaved 
method. This is indeed the case. Figure [7] also shows re¬ 
sults with coupled cl uster w ith single, double, and triple 
excitations (CCSDT) j 43 l 44 l as well as with what we call 
CCSDTO, by which we mean CCSDT subject to the re¬ 
striction that T 2 be singlet paired as in CCDO. Clearly, 
CCSDT has the same qualitative failures as do CCD 
and CCSD, but once the appropriate modifications are 
made to the T 2 operator, the resulting method is again 
well-behaved. Adding correlations atop CCDO, in other 
words, is fairly straightforward. The only challenge is ap¬ 
propriately correcting for the triplet-paired portion T .^ 
of T 2 which CCDO has omitted; this portion is clearly 
important for getting accurate results in regions where 
restricted CCD works well, but must be included care¬ 
fully so as not to reproduce the failures of CCD as well 
as its successes. 
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Rh-h(A) 



FIG. 8 . Total energies for the H 14 ring as a function of inter¬ 
atomic separation. The curves marked “SA-RHF” and “SB- 
RHF” are RHF solutions which respectively preserve or break 
spatial symmetry, while “SA-CCSDO” and “SB-CCSDO” are 
CCSDO curves based on the SA-RHF and SB-RHF references. 
The SA-CCSD curve is CCSD based on the SA-RHF refer¬ 
ence. We could not converge CCSD based upon the SB-RHF 
determinant. The black points show the onset of instability of 
the SA-RHF determinant toward UHF or the SB-RHF state. 


A. Instabilities in CCDO 

Thus far, we have shown no examples for which CCDO 
breaks down. In fact, while CCDO is much less prone to 
failure for strongly correlated systems than is CCD, it is 
not impervious to these problems. When problems are 
sufficiently strongly correlated, that is to say, CCDO is 
a great improvement upon CCD but is not necessarily a 
panacea. 

Figure [8] shows results for a ring of fourteen equally- 
spaced hydrogen atoms in the STO-3G basis set. As 
the radius of the ring increases, the spacing between ad¬ 
jacent atoms increases as well. Loosely, this forms an 
analogy with the Hubbard Hamiltonian, though in the 
Hubbard Hamiltonian the electron-electron interaction is 
local and has zero range. Even before equilibrium, there 
is an instability in the symmetry-adapted RHF (marked 
as “SA-RHF” on the plot) toward a UHF solution. For 
larger bond lengths, the symmetry-adapted RHF devel¬ 
ops an additional instability, toward another RHF so¬ 
lution which breaks spatial symmetry (marked as “SB- 
RHF”) on the plot. This symmetry-broken RHF has a 
dinrer-like structure, and the occupied orbitals essentially 
resemble seven sets of H 2 bonding orbitals. 

Conventional CCSD based on the symmetry-adapted 
RHF reference begins to overcorrelate dramatically 
shortly past equilibrium, after which we can no longer 
converge the CCSD equations. We were unable to 
converge the CCSD equations on the symmetry-broken 
RHF reference. In contrast, we were able to converge 


FIG. 9. Total energies in the 10-site Hubbard Hamiltonian at 
half filling as a function of interaction strength U/t. 


CCSDO based on the symmetry-adapted reference until 
much larger bond lengths. At the point for which the 
symmetry-adapted RHF develops an instability toward 
spatial symmetry breaking, the CCSDO begins to break 
down. We can, however, solve this problem by using the 
symmetry-broken RHF as a reference for CCSDO. The 
same general observations hold for CCD and CCDO. 

Because this hydrogen ring is loosely analogous to 
the Hubbard Hamiltonian, it may not be too surprising 
that the Hubbard Hamiltonian displays a similar phe¬ 
nomenon. In Fig. [9] we show results for the 10-site 
Hubbard Hamiltonian. We see that CCD fails dramati¬ 
cally and CCDO begins to break down for large enough 
U/t. This is, in fact, also true for the smaller 6-site Hub¬ 
bard problem, where the CCD energy reaches a maxi¬ 
mum near U/t ss 6.2 and the CCDO energy is maximal 
for U/t ss 20; eliminating the ring and crossed-ring con¬ 
tractions entirely, in contrast, undercorrelates badly but 
leaves a total energy which increases monotonically at 
least until U/t = 30. We have not yet been able to ex¬ 
plain the origin of the failure of CCDO for this problem 
although such large U values are normally outside the 
regime of physical systems. The RHF solution we have 
used as a reference is stable toward spatial symmetry 
breaking, unlike with the hydrogen ring, so it is not clear 
how to use a different RHF reference to resolve this prob¬ 
lem. We have traced the failure of CCDO in this case to 
the ring diagrams, as removing the ring and crossed-ring 
terms leaves, as one might expect, a well-behaved but 
inaccurate result. 


B. Triplet-Paired CCD 

Thus far, we have shown that by simply choosing the 
pair creation and annihilation operators defining the T 2 
cluster operator to be singlet paired only, we obtain a 
method which is largely free of the major failures of 
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FIG. 10. Total energies of the N 2 molecule in the STO-3G 
basis set as a function of bond length. 


CCD. One could alternatively have chosen the pair cre¬ 
ation and annihilation operators to be triplet paired. In 
other words, instead of retaining only T .^ we could have 
retained instead only T^\ This defines what we call 
CCD1, where the “1” indicates the 5 = 1 triplet-pairing 
nature of the correlator. Where CCDO entirely neglects 
correlations between same-spin electrons, such correla¬ 
tions are present in CCD1, though at the price of less 
focus on the more dominant opposite-spin correlations. 
We emphasize again that the classification of correlations 
into “same spin” and “opposite spin” is different from 
the classification into singlet- and triplet-paired corre¬ 
lations. Singlet pairing is opposite-spin only but does 
not contain all opposite-spin correlations; the remain¬ 
der of the opposite-spin correlations, together with all 
same-spin correlations, are contained within the triplet 
pairing channel. In Fig. [l0]we show the dissociation of 
N 2 in the minimal basis one more time, now including 
CCD1. As one can see, CCD1, like CCDO, is generally 
stable and undercorrelating. This is also clearly true in 
the Hubbard Hamiltonian, for which CCD1 predicts zero 
correlation energy (because the two-electron interaction 
in the Hubbard Hamiltonian is singlet paired only). We 
emphasize that RHF-based CCD has a triplet pairing 
channel for which the three distinct components all have 
the same wave function amplitudes. In UHF-based CCD, 
one could attempt to define a triplet pairing channel but 
due to broken spin symmetry would obtain three distinct 
sets of T 2 amplitudes. 

Individually, CCDO and CCD1 are generally well- 
behaved and undercorrelating. It is apparent that one 
cannot meaningfully just solve the CCDO equations for 
the symmetric part of the cluster amplitude and 
the CCD1 equations for the antisymmetric part; in the 
Hubbard Hamiltonian, this would return the CCDO re¬ 
sult, while in N 2 it would overcorrelate wildly. The cou¬ 
pling between these singlet- and triplet-paired parts of T 2 


which is present in the full CCD equations thus causes 
CCD to deliver reasonable energies in weakly correlated 
systems; however, the interaction between singlet- and 
triplet-paired channels also apparently causes the failures 
of CCD for strongly correlated systems. It seems at least 
possible, then, that modifying the coupling between these 
two parts of T 2 might provide a useful alternative, albeit 
one unexplored here. It is evident that full inclusion of 
higher cluster amplitudes (T 3 , T 4 , and so on) serves to 
renormalize the coupling between the T 2 amplitudes. 


IV. CONCLUSIONS 

The failure of traditional coupled cluster theory in 
strongly correlated systems is well known, and is usu¬ 
ally understood as a failure of the reference determinant 
to even qualitatively describe the wave function. When 
such failures take place, they can be ameliorated in one 
of two ways: either the reference must be improved, or 
the correlator used to create the correlated wave function 
from the reference must be improved. Neither option, un¬ 
fortunately, is entirely ideal. Improving the reference can 
be done by adopting some form of multireference coupled 
cluster approach, but this is by no means straightforward. 
Improving the correlator by making it more complete - 
that it, by approaching nearer to full configuration in¬ 
teraction - is certainly conceptually straightforward, but 
computationally can be very demanding. Typically we 
perforce resort to the use of a broken symmetry mean- 
field reference; for repulsive two-body interactions, spin 
symmetry is broken, while for attractive two-body inter¬ 
actions, we break number symmetry instead.® 

What has been explored less is improving the correla¬ 
tor by removing pieces. This has the advantage of com¬ 
putational simplicity but the disadvantage that what is 
to be removed must be chosen with a certain amount 
of care. The point of this manuscript is simply that in 
CCDO we have found one way to simplify the cluster op¬ 
erator so as to at least substantially remove the failures 
of the cluster operator of CCD in the strongly correlated 
limit. By working directly at the level of modifying the 
cluster operator, CCDO ensures that we can follow all 
the usual machinery of standard coupled cluster theory 
straightforwardly. In particular, we have a well-defined 
fermionic wave function. Just as with standard CCD, we 
can define a left-hand wave function (0|(1 + Z) exp(—T) 
where Z is a de-excitation operator which has, in CCDO, 
symmetric coefficients z^ b and with the aid of which we 
can construct density matrices and take expectation val¬ 
ues. We can in principle define excited states via equation 
of motion coupled cluster, and so on. 

Of course the pieces of the cluster operator that CCDO 
removes must be put back in somehow. That is, we do not 
advocate CCDO as a stand-alone method for the descrip¬ 
tion of electronic systems, because the same-spin corre¬ 
lations have been removed from the theory entirely and 
results in the absence of strong correlation are poorer 
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than what we could get from CCD itself. Moreover, the 
opposite-spin correlation in CCDO is not the same as in 
CCD either, as can be seen in the Hubbard Hamiltonian 
where there is no same-spin interaction, yet CCD and 
CCDO deliver different results nonetheless. One must, at 
the end of the day, correct CCDO for what it has taken 
out of CCD. 

Let this not, however, obscure the point of this paper: 
to the extent that one can identify which terms in the 
coupled cluster equations are responsible for the failures 
of the method, the failures can be resolved cheaply by 
simply removing those terms in a suitable way, or more 
expensively by adding additional terms which counteract 
the effect. The former approach is computationally sim¬ 
pler and, we feel, may merit more attention than it has 
thus far received. Identifying the culprit at the CCD level 
may help us design a selective introduction of higher clus¬ 
ter operators to balance their defects. Work along these 
lines will be presented in due time. 

By virtue of its exponential parameterization, single¬ 
reference coupled cluster theory can in principle deliver 
the correct coefficients on several important determinants 
without breaking down entirely, if we choose the correla¬ 


tor correctly; in other words, while based on a single ref¬ 
erence determinant, the wave function can be multicon- 
figurational in nature. It does not possess full flexibility 
to choose the coefficients on each determinant indepen¬ 
dently, but this is apparently not necessary to describe 
static correlation in many situations. The successes of 
coupled cluster methods which describe strongly corre¬ 
lated systems by removing pieces of the cluster operator 
provide a fresh perspective on static correlation. Cou¬ 
pled cluster methods which are based on a symmetry- 
adapted reference which is unstable toward broken sym¬ 
metry can at least partially account for static correlation 
effects whenever the coupled-cluster energy is stable and 
well-behaved. 
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